Assignment #5

# SEE modeldata package for new datasets
library(tidyverse)         # for graphing and data cleaning
library(tidymodels)        # for modeling
library(ranger)            # for random forest - will need for shiny app
library(lubridate)         # for date manipulation
library(themis)            # for up and downsampling
library(DALEX)             # for model interpretation  
library(DALEXtra)          # for extension of DALEX
library(lime)
theme_set(theme_minimal()) # Lisa's favorite theme

data("lending_club")
# Data dictionary (as close as I could find): https://www.kaggle.com/wordsforthewise/lending-club/discussion/170691

Put it on GitHub! – DONE

Github Repo Link: https://github.com/apalma127/assignment-5

Interpretable ML methods – DONE

We will once again use the lending club data that we used in the 3rd assignment. We will focus on the random forest model, which I recreate below. (Note we use this model even though the true negative rate of the training set is quite bad.)

set.seed(494) # for reproducibility

#split data
lending_split <- initial_split(lending_club,
                               prop = .75,
                               strata = Class)

lending_training <- training(lending_split)
lending_test <- testing(lending_split)


#create recipe - including up and downsampling for model fitting
set.seed(456)
rf_recipe <- 
  recipe(Class ~ .,
         data = lending_training) %>% 
  step_upsample(Class, over_ratio = .5) %>% 
  step_downsample(Class, under_ratio = 1) %>% 
  step_mutate_at(all_numeric(), 
                 fn = ~as.numeric(.))

# create model
rf_model <- 
  rand_forest(mtry = tune(), 
              min_n = tune(), 
              trees = 100) %>% 
  set_mode("classification") %>% 
  set_engine("ranger")

# create workflow
rf_workflow <-
  workflow() %>% 
  add_recipe(rf_recipe) %>% 
  add_model(rf_model)

  grid_regular(finalize(mtry(),
                        lending_training %>%
                          select(-Class)),
               min_n(),
               levels = 3)

# create penalty grid
  rf_penalty_grid <- 
grid_regular(finalize(mtry(),
                        lending_training %>%
                          select(-Class)),
               min_n(),
               levels = 3)


# create cv samples
set.seed(494) #for reproducible 5-fold
lending_cv <- vfold_cv(lending_training,
                       v = 5)

# tune model
rf_tune <- 
  rf_workflow %>% 
  tune_grid(
    resamples = lending_cv,
    grid = rf_penalty_grid
  )

# find model with best accuracy
best_accuracy <-
  rf_tune %>% 
  select_best(metric = "accuracy")

# finalize model
rf_final <- rf_workflow %>% 
  finalize_workflow(best_accuracy) %>% 
  fit(data = lending_training)

1. Use functions from the DALEX and DALEXtra libraries to create a histogram and boxplot of the residuals from the training data.

rf_explain <- 
  explain_tidymodels(
    model = rf_final,
    data = lending_training %>% select(-Class), 
    y = lending_training %>% 
      mutate(Class_num = as.integer(Class =="good")) %>% 
      pull(Class_num),
    label = "rf"
  )

## Preparation of a new explainer is initiated
##   -> model label       :  rf 
##   -> data              :  7392  rows  22  cols 
##   -> data              :  tibble converted into a data.frame 
##   -> target variable   :  7392  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.3 , task classification (  default  ) 
##   -> predicted values  :  numerical, min =  0 , mean =  0.8477854 , max =  1  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.095 , mean =  0.1004018 , max =  0.72  
##   A new explainer has been created! 

rf_mod_perf <-  model_performance(rf_explain)

rf_mod_perf

## Measures for:  classification
## recall     : 0.9897275 
## precision  : 1 
## f1         : 0.9948372 
## accuracy   : 0.9902597 
## auc        : 1
## 
## Residuals:
##     0%    10%    20%    30%    40%    50%    60%    70%    80%    90%   100% 
## -0.095  0.000  0.010  0.020  0.040  0.060  0.080  0.120  0.180  0.270  0.720

hist_plot <- 
  plot(rf_mod_perf,
       geom = "histogram")


box_plot <-
  plot(rf_mod_perf,
       geom = "boxplot")

hist_plot

box_plot

How do they look? Any interesting behavior?

They are not centered around 0 as you would want a good model to do. They appear to skew right heavily well into the positives indicating a continuous underprediction. This is very significant as it shows a pattern of underpredictions which shouldn’t be a constant pattern if it is a good model.

2. Use DALEX functions to create a variable importance plot from this model.

rf_var_imp <- 
  model_parts(
    rf_explain
    )

plot(rf_var_imp, show_boxplots = TRUE)

What are the most important variables?

The most importannt variables used to help correctly predict if a loan was paid back was interest rate, sub_grade, open_il_24m, and annual income. These are not surprising at all because most of them relate exactly back to affecting one’s ability to repay … ie if an interest rate is incredibly high it makes the pay back extremely less likely, if one has a very low income it makes payback very difficult….

3. Write a function called cp_profile to make a CP profile.

The function will take an explainer, a new observation, and a variable name as its arguments and

create a CP profile for a quantitative predictor variable.

You will need to use the predict_profile() function inside the function you create - put the variable name there so the plotting part is easier.

obs2 <- lending_training %>% 
  slice(2)
obs2

CODE TEST FOR ANNUAL INC AND INT RATE

is.integer(lending_training$annual_inc)

## [1] FALSE

is.integer(lending_training$int_rate)

## [1] FALSE

TEST 1

cp_profile_test <- predict_profile(explainer = rf_explain, 
                          new_observation = obs2,
                          variables = c("annual_inc", "int_rate"))


cp_profile_test

cp_profile_test %>% 
  filter(`_vname_` %in% c("annual_inc")) %>% 
  ggplot(aes(x = annual_inc,
             y = `_yhat_`)) +
  geom_line() 

With y_hat being our predicted class outcome (with 1 = good pay back on time and 0 = bad not paid back on time) it appears to show us: as income increases, while it doesn’t appear to be a major difference, it is slightly more likely the person is not on time w repayment. Yet, it appears it is much more constant w the trend line for higher income whereas lower income is much more chaotic in predicted payback.

TEST 2

cp_profile_test %>% 
  filter(`_vname_` %in% c("int_rate")) %>% 
  ggplot(aes(x = int_rate,
             y = `_yhat_`)) +
  geom_line() 

With y_hat being our predicted class outcome (with 1 = good pay back on time and 0 = bad not paid back on time) it appears to show us a very obvious trend. As expected, as the interest rate increases meaning increasingly more money is owed back than borrowed, the ability to pay on time and the predicted class value races downward closer to bad meaning the payment is much less likely to be on time.

Now Function Writing time…

Write a function called cp_profile to make a CP profile.

The function will take an explainer, a new observation, and a variable name as its arguments and

create a CP profile for a quantitative predictor variable.

You will need to use the predict_profile() function inside the function you create - put the variable name there so the plotting part is easier.

You’ll also want to use .data[[]] rather than aes() and quote the variables. Use the cp_profile() function to create one CP profile of your choosing.

Be sure to choose a variable that is numeric, not integer. There seem to be issues with those that I’m looking into.

 cp_profile <- function(explainer_fn, new_obs, `var`) {
  
  
  profile <- predict_profile(explainer = explainer_fn, 
                          new_observation = new_obs,
                          variables = `var`)
  
  graph <- profile %>% 
    filter(`_vname_` %in% c(`var`)) %>% 
    ggplot(aes(x = .data[[`var`]], y = `_yhat_`)) +
    geom_line() 

  
  graph
 }

lending_numeric <- select_if(lending_training, is.numeric)            

lending_numeric

cp_profile(rf_explain, obs2, "annual_inc")

cp_profile(rf_explain, obs2, "int_rate")

cp_profile(rf_explain, obs2, "funded_amnt")

cp_profile(rf_explain, obs2, "revol_util")

4. Use DALEX functions to create partial dependence plots (with the CP profiles in gray) for the 3-4 most important variables.

If the important variables are categorical, you can instead make a CP profile for 3 observations in the dataset and discuss how you could go about constructing a partial dependence plot for a categorical variable (you don’t have to code it, but you can if you want an extra challenge).

Most Important: interest rate, sub_grade, open_il_24m, and annual income

is.factor(lending_training$int_rate)

## [1] FALSE

is.factor(lending_training$sub_grade)

## [1] TRUE

is.factor(lending_training$open_il_24m)

## [1] FALSE

is.factor(lending_training$annual_inc)

## [1] FALSE

rf_pdp <- model_profile(explainer = rf_explain, 
                        variables = c("int_rate", "open_il_24m", "annual_inc"))

plot(rf_pdp, 
     variables = c("int_rate", "open_il_24m", "annual_inc"),
     geom = "profiles")

Discuss how you could go about constructing a partial dependence plot for a categorical variable

– Subgrade = Categorical

– Steps to Create PDP for subgrade

++++ one major step to do and then can use code above: convert from categorical / factor to numeric value !!!

– for each level of subgrade, working from top to bottom, assign a numeric value using a series of ifelse statements within mutate for subgrade within lending training

– now, w subgrade being numeric, can feed right into the list of variables above…

5. Choose 3 observations and do the following for each observation:

OBSERVATION 2

obs2_q5 <- lending_test %>% 
  slice(2)
obs2_q5

• Construct a break-down plot using the default ordering.

pp_rf_2 <- predict_parts(explainer = rf_explain,
                          new_observation = obs2_q5,
                          type = "break_down") 

plot(pp_rf_2)

Interpret the resulting graph. Which variables contribute most to each observation’s prediction?

First we can seethe initial bars start at the intercept of 0.848 which marks the value of the average predicted class (1 = good, 0 = bad) when applying the rf model to the training data.

We can see we end up 0.108 lower in our predicted class of paying back a loan than the intercept thanks to some heavy hitting negative variables. The resulting 0.74 is a fairly strong good prediction.

We can see that the largest and most influential tugs from variables down towards the actual predicted value from the intercept comes from:

total_bal_il == 55445 –> -0.048

num_il_tl == 10 –> -0.032

open_il_24m == 2 –> -0.029

This means that for each of these variables, if they were fixed at the values they are set to, the change in average class prediction would be these large negative values.

• Construct a SHAP graph and interpret it.

rf_shap_2 <-predict_parts(explainer = rf_explain,
                        new_observation = obs2_q5,
                        type = "shap",
                        B = 20 
)

plot(rf_shap_2)

Does it tell a similar story to the break-down plot?

We can see that SEVERAL variables have boxplots of their effects straddling zero spreading across both positive and negative values.

These variables set to their appropriate constants include:

total_bal_il

int_rate

sub_grade

num_il_tl

emp_length

total_il_high_credit_limit

This means that when we changed the order of considering these variables and re-run the breakdown test 20 times, we have the variables having multiple instances of both negative and positive effects.

This narrative is in direct conflict with the clear cut picture we were painted above about variable effects.

• Construct a LIME graph (follow my code carefully).

model_type.dalex_explainer <- DALEXtra::model_type.dalex_explainer
predict_model.dalex_explainer <- DALEXtra::predict_model.dalex_explainer

set.seed(2)
lime_rf_2 <- predict_surrogate(explainer = rf_explain,
                             new_observation = obs2_q5 %>%
                               select(-Class), 
                             n_features = 5,
                             n_permutations = 1000,
                             type = "lime")

lime_rf_2 %>% 
  select(model_r2, model_prediction, prediction) %>% 
  distinct()

plot(lime_rf_2) +
  labs(x = "Variable")

How close is each original prediction to the prediction from the local model? Interpret the result. You can also try using fewer or more variables in the local model than I used in the example.

Prediction from Original RF Model: 0.74

Prediction from Local Model: 0.82 (+0.08)

The predictions are pretty close…the Local Model is +0.08…

– They are basically the same because most variables have roughly the same impact… probably because the total_il_high_credit_limit variable appears to have slightly greater impact here and some variables impacts may be on the higher end of their boxplots…also just so happens 3 / 5 vars have positive impacts of all same size….

Local Model R^2: 0.05 – Not Good

Bar Plot Shows Variable Importance From Local Model:

– annual income nd total bal il have big negative effects

– the rest appear to have equally as high impact vars but positive

OBSERVATION 20

obs20 <- lending_test %>% 
  slice(20)
obs20

• Construct a break-down plot using the default ordering.

pp_rf_20 <- predict_parts(explainer = rf_explain,
                          new_observation = obs20,
                          type = "break_down") 

plot(pp_rf_20)

Interpret the resulting graph. Which variables contribute most to each observation’s prediction?

First we can see the initial bars start at the intercept of 0.848 which marks the value of the average predicted class (1 = good, 0 = bad) when applying the rf model to the training data.

We can see we end up 0.018 lower in our predicted class of paying back a loan than the intercept. This very little change is due to the fact that a lot of the variables’ effects have cancelling effects (fairly even split of + and - effects). The resulting 0.83 is a strong good prediction.

We can see that the largest tugs come from variables both + and - (and I guess that makes them most influential as they cancel each other out for the most part)…:

int_rate == 12.99 –> -0.04

total_bal_il == 26275 –> 0.037

inq_fi == 4 –> -0.028

This means that for each of these variables, if they were fixed at the values they are set to, the change in average class prediction would be their outputted vals.

This was a weird one to eval contribution: 1) do the bigger tugs still count as most significant because they are cancelling eachother out or 2) do the small negative tugs count more even though they wouldn’t matter if the big tugs didn’t cancel out????

• Construct a SHAP graph and interpret it.

rf_shap_20 <-predict_parts(explainer = rf_explain,
                        new_observation = obs20,
                        type = "shap",
                        B = 20 
)

plot(rf_shap_20)

Does it tell a similar story to the break-down plot?

We can see that SEVERAL variables have boxplots of their effects straddling zero spreading across both positive and negative values.

These variables set to their appropriate constants include:

total_bal_il

sub_grade

annual_inc

int_rate

addr_state

all_util

inq_last_6mths

This means that when we changed the order of considering these variables and re-run the breakdown test 20 times, we have almost every single variable having multiple instances of both negative and positive effects.

This narrative is in direct conflict with the clear cut picture we were painted above as the vars effects are not as clear as we thought above.

• Construct a LIME graph (follow my code carefully).

model_type.dalex_explainer <- DALEXtra::model_type.dalex_explainer
predict_model.dalex_explainer <- DALEXtra::predict_model.dalex_explainer

set.seed(2)
lime_rf_20 <- predict_surrogate(explainer = rf_explain,
                             new_observation = obs20 %>%
                               select(-Class), 
                             n_features = 5,
                             n_permutations = 1000,
                             type = "lime")

lime_rf_20 %>% 
  select(model_r2, model_prediction, prediction) %>% 
  distinct()

plot(lime_rf_20) +
  labs(x = "Variable")

How close is each original prediction to the prediction from the local model? Interpret the result. You can also try using fewer or more variables in the local model than I used in the example.

Prediction from Original RF Model: 0.83

Prediction from Local Model: 0.8 (-0.03)

The predictions are fairly close…the Local Model is -0.03…

– This likely is because of the massive negative impact of revol_util and annual_inc which pushed the pred slightly below as they countered the pos impact from others. Annual Income’s effects appear to be much more negative here than above…

Local Model R^2: 0.07 – Terrible

Bar Plot Shows Variable Importance From Local Model:

– revol_util has a massive negative impact

– annual income has equally massive negative impact

– inq last 6 mths and total il high credit have moderately large positive impacts

OBSERVATION 200

obs200 <- lending_test %>% 
  slice(200)
obs200

• Construct a break-down plot using the default ordering.

pp_rf_200 <- predict_parts(explainer = rf_explain,
                          new_observation = obs200,
                          type = "break_down") 

plot(pp_rf_200)

Interpret the resulting graph. Which variables contribute most to each observation’s prediction?

First we can see the initial bars start at the intercept of 0.848 which marks the value of the average predicted class (1 = good, 0 = bad) when applying the rf model to the training data.

We can see we end up 0.132 HIGHER in our predicted class of paying back a loan than the intercept thanks to some big positive jumps from influential variables. The resulting 0.98 is an extremely strong good prediction.

We can see that the largest positive jumps from variables come from:

int_rate == 13.44 –> 0.029

sub_grade == 13 –> 0.029

total_bal_il == 0 –> 0.028

and also it helps that all but 3 in this case have positive impacts…

This means that for each of these variables, if they were fixed at the values they are set to, the change in average class prediction would be these large positive values which pulled our prediction for this observation up very high.

• Construct a SHAP graph and interpret it.

rf_shap_200 <-predict_parts(explainer = rf_explain,
                        new_observation = obs200,
                        type = "shap",
                        B = 20 
)

plot(rf_shap_200)

Does it tell a similar story to the break-down plot?

We can see that SEVERAL variables have boxplots of their effects straddling zero spreading across both positive and negative values.

These variables set to their appropriate constants include:

total_bal_il

sub_grade

int_rate

total_il_high_credit_limit

annual_inc

num_il_til

This means that when we changed the order of considering these variables and re-run the breakdown test 20 times, we have these variables having multiple instances of both negative and positive effects.

This narrative is in direct conflict with the clear cut picture we were painted above as we see the impacts of vars are not as obvious.

• Construct a LIME graph (follow my code carefully).

model_type.dalex_explainer <- DALEXtra::model_type.dalex_explainer
predict_model.dalex_explainer <- DALEXtra::predict_model.dalex_explainer

set.seed(2)
lime_rf_200 <- predict_surrogate(explainer = rf_explain,
                             new_observation = obs200 %>%
                               select(-Class), 
                             n_features = 5,
                             n_permutations = 1000,
                             type = "lime")

lime_rf_200 %>% 
  select(model_r2, model_prediction, prediction) %>% 
  distinct()

plot(lime_rf_200) +
  labs(x = "Variable")

How close is each original prediction to the prediction from the local model? Interpret the result. You can also try using fewer or more variables in the local model than I used in the example.

Prediction from Original RF Model: 0.98

Prediction from Local Model: 0.81 (-0.17)

The predictions are not close…the Local Model is -0.17…

– This likely is because here we have massive negative impacts seen from revol_util and annual_income, much higher than displayed above, which account for the massive drop in predicted class score.

Local Model R^2: 0.06 – Awful

Bar Plot Shows Variable Importance From Local Model:

– revol_util, annual income have massively large negative impacts that outpace the positive impacts of inq last 6mths, inq_fi, and total_bal_il.

6. Describe how you would use the interpretable machine learning tools we’ve learned (both local and global) in future machine learning projects? How does each of them help you?

I could find great use in using the boxplot and histogram of residuals for a model from global interpretable ML and from local interpretable ML I would definitely have great value from using breakdown profiles and shapley plots.

For Global: Visualizing Residuals

This past Summer, I had an internship with a Fintech company called AvidXChange, a Charlotte, NC based Accounts Payable Automation company. I used sql to pull data from databases and machine learning skills to use neural networks and deep learning. We wanted to more accurately predicted which of our clients accepted Virtual Credit Card and which didn’t as this was a big problem. Avid makes most of their money from fees on VCC cards being used to pay off invoices as opposed to just direct deposit (this is more ideal for most as they get the money much quicker despite a very small fee). If able to better identify who simply can’t take VCC, it would save a lot of time and money better used on other clients.

When we finalized the mode, I would have loved to see how well it predicted through seeing the distribution of its residuals. My model also was a class prediction like w the lending data, specifically using the spelling and make up of the company name to predict if they took VCC (the hunch was that companies that appeared to be people did not accept VCC and the model used a series of dictionaries for letters and names from census data to distinguish names).

If I could analyze the residuals a bit more with things like the boxplot and histogram like we did above, I would have been able to see a bit better how well we predicted correctly. Specially, I would have been able to see if it was a good model beyond simply looking at output metrics like I did. If I could have visualized the residuals and seen how it tended to predict the majority class more than most other things, that would have been very easy to see and cool to help show the results as opposed to simply telling people the results with metrics filled with jargon.

For Local: Breakdown Profiles and Shapley

Also, in my internship I added on to a VCC model a propensity score value which predicted a probability at which a company (we thought) would accept VCC. It took in many variables including past invoices, size of the invoice, size of the company, if they have taken VCC in the past before, location…

I would have loved to see for specific companies which variables were more significant in driving the propensity scores to where they went. Obviously Panera takes VCC because they take credit card but what variable specifically was honing in on this?

Seeing the little tug of war between variables and how they effected the outcome woud have been really cool. Additionally, I would have liked to see the permanence of these variables effects beyond just one observation. If I could have also used a shap plot, I could have seen how firm these variable impacts were as well and seen how it changed from large to small company occurrences.

7. Save this final model using the write_rds() function - see the section of the tidymodels intro for a similar example, but we’re using write_rds() instead of saveRDS(). We are going to use the model in the next part. You’ll want to save it in the folder where you create your shiny app. Run the code, and then add eval=FALSE to the code chunk options (next to the r inside the curly brackets) so it doesn’t rerun this each time you knit.

# # finalize model
# rf_final <- rf_workflow %>% 
#   finalize_workflow(best_accuracy) %>% 
#   fit(data = lending_training)


write_rds(rf_final, "rf_final.rds")

rf_final_read <- readRDS("rf_final.rds")

Shiny app – DONE

My App: https://ajp28117.shinyapps.io/Lending_Club/

Repo for App: https://github.com/apalma127/lending_small

Website Post: https://anthonypalma.netlify.app/posts/shinyapps/ (~~ scroll down to lending data shiny app ~~)

Data Ethics: Data visualization principles – DONE

Were there any principles mentioned that you hadn’t heard of before?

I was unaware that the y axis on line graphs shouldn’t include 0 and that the graph should be much more honed in on the region of interest than zoomed out. I actually have never done this to be honest just because I thought it was good etiquette to always start at 0…

What graph stood out for you as “the worst”?

The worst has to be the bubble chart with the professor in front. While it does have a lot of data and cool things going on, it is a mess. With no axes and random size differences, there is so much going on the only thing you can see is relative difference. Additionally, there is a lot of blobbing of points on top of each other making even relative difference difficult to see.

Did any of the graphs fool you?

The gun deaths in FL graph was veryyyyy confusing because of the upside down filling of it. At first I was like wow super weird stand your ground DECREASED gun deaths but that was simply because the creator of the graph wanted you to see quite the opposite of reality.

How does practicing good data visualization principles fit in with data ethics?

The florida gun deaths one was a really good example of how bad data practices, no matter the intent, can actually manipulate the truth and lie to the viewer. Such fake news would be damaging to a narrative in an election especially because of how contrary to reality it is. If someone were to do this during election time with important stats about abortion and immigration, it would be crucial to changing minds and votes.